%--------------------------------------------------------------------------
% CDS PRICING
%--------------------------------------------------------------------------


% =========================================================================
% default probabilities under P
% =========================================================================

[Cheb_nodes,~] = qnwcheb(nE,-1,1);
k_vec       = repmat(k',[nS,1]);
k_vec       = k_vec(:); 
kb_vec      = repmat(k_bar,[nK,1]);
kb_vec      = kb_vec(:); % [nS*nK,1]

mu_vec      = repmat(mu_X,[nK,1]);
mu_vec      = mu_vec(:); % [nS*nK,1]
sig_vec     = repmat(sig_X,[nK,1]);
sig_vec     = sig_vec(:); % [nS*nK,1]
P_mat       = repmat(P,[nK,1]); % [nS*nK,nS]

default     = reshape( k' <= k_def ,[nK*nS,1]);
issuance    = reshape( k' >= k_iss ,[nK*nS,1]);

Amat        = (log(k_def)' - mu_vec - log(k_vec) )./sig_vec;
Amat_iss    = (log(k_def)' - mu_vec - log(kb_vec))./sig_vec;
Amat(issuance,:) = Amat_iss(issuance,:);

p1          = reshape(~default.*sum( normcdf(Amat) .* P_mat ,2),[nS,nK]); % [nS,nK]
p           = zeros(nS,nK,Tmax_CDS);
p(:,:,1)    = p1;

for t=2:Tmax_CDS % maturity
    for is=1:nS % future state
        a       = max(-5,(log(k_def(is)./k_vec) - mu_vec)./sig_vec); 
        b       = 5;     
        Eshocks = (b-a)*(Cheb_nodes'+1)/2 + a;
        temp    = pi*(b-a)/(2*nE);
        weights = temp.*sqrt(1-Cheb_nodes'.^2);
        probE   = normpdf(Eshocks).*weights;
        
        % kappa transition
        kp             = k_vec  .* exp(mu_vec + sig_vec.*Eshocks); %
        kp_issuance    = kb_vec .* exp(mu_vec + sig_vec.*Eshocks); %
        kp(issuance,:) = kp_issuance(issuance,:);
        kp(default,:)  = k_min; % any default prob evaluated at k_min equals zero
        
        p_Fit    = griddedInterpolant(k,p(is,:,t-1)','linear');
        p(:,:,t) = p(:,:,t) + reshape(P_mat(:,is) .* sum(p_Fit(kp).*probE,2),[nS,nK]); % [nS,nK] %%%
    end
end


% =========================================================================
% loss given default under P
% =========================================================================


lamD_fit        = griddedInterpolant({1:nS,k},lamD,'linear');

mu_vec          = repmat(mu_X,[nK*nS,1]);
sig_vec         = repmat(sig_X,[nK*nS,1]);

k_vec           = repmat(k',[nS,1,nS]);
k_vec           = k_vec(:);

VEC_k_bar       = reshape(repmat(k_bar',[nK*nS,1]),[nK*nS^2,1]);
kb_vec          = repmat(k_bar',[nK*nS,1]);
kb_vec          = kb_vec(:);

VEC_lamD        = reshape(repmat(lamD_fit((1:nS)',k_bar)',[nK*nS,1]),[nK*nS^2,1]);

MAT_P           = reshape(permute(repmat(P,[nK,1,nS]),[1,3,2]),[nK*nS^2,nS]);
default         = reshape(repmat(k'<=k_def,[1,1,nS]),[nK*nS^2,1]);
issuance        = reshape(repmat(k'>=k_iss,[1,1,nS]),[nK*nS^2,1]);

Amat            = (log(k_def)' - mu_vec - log(k_vec) )./sig_vec;
Amat_iss        = (log(k_def)' - mu_vec - log(kb_vec))./sig_vec;
Amat(issuance,:) = Amat_iss(issuance,:);

MAT_cdf         = normcdf( Amat-sig_vec );
MAT_Rf          = permute(repmat(Rf(:,1:Tmax_CDS),[1,1,nK,nS]),[1,3,4,2]);
L1              = repmat(p1(:),[nS,1]) - MAT_P .* MAT_cdf .* k_vec .* exp(mu_vec+sig_vec.^2/2) ./ (VEC_lamD.*VEC_k_bar) * ((1-omega).*lamA);
L1              = max(L1,0);
L1              = reshape(L1,[nS,nK,nS]);

% transition of most recent issuance state (note: predetermined)
markov_state     = repmat((1:nS)',[nK*nS,1]);
xi_s_p           = reshape(repmat(1:nS,[nK*nS,1]),[nK*nS^2,1]);
xi_s_p(issuance) = markov_state(issuance);
xi_s_p_dummy     = zeros(nK*nS^2,nS);
xi_s_p_dummy((xi_s_p-1)*nK*nS^2 + (1:nK*nS^2)') = 1;

% multi-period L
L_p             = zeros(nS,nK,nS,Tmax_CDS);
L_p(:,:,:,1)    = L1;

for t=2:Tmax_CDS % maturity
    
    L_temp = zeros(nS^2*nK,nS); % 2nd dim: future last issuace state
    
    for is=1:nS % future markov state
        
        a       = max(-5,(log(k_def(is)./k_vec) - mu_vec)./sig_vec); 
        b       = 5;        
        Eshocks = (b-a)*(Cheb_nodes'+1)/2 + a;
        temp    = pi*(b-a)/(2*nE);
        weights = temp.*sqrt(1-Cheb_nodes'.^2);
        probE   = normpdf(Eshocks).*weights;
        
        % kappa transition
        kp          = k_vec  .* exp(mu_vec + sig_vec.*Eshocks); %
        kp_issuance = kb_vec .* exp(mu_vec + sig_vec.*Eshocks); %
        kp(issuance,:) = kp_issuance(issuance,:);
        kp(default,:)  = k_min; % any default prob evaluated at k_min equals zero  
        
        for is_s=1:nS % last issuance state
            L_Fit          = griddedInterpolant(k,squeeze(L_p(is,:,is_s,t-1)),'linear');
            L_temp(:,is_s) = L_temp(:,is_s) + MAT_P(:,is) .* nansum(L_Fit(kp).*probE,2);
        end
    end
    L_p(:,:,:,t) = reshape(nansum(L_temp.*xi_s_p_dummy,2),[nS,nK,nS]);
end

LGD_p   = L_p ./ permute(repmat(p,[1,1,1,nS]),[1,2,4,3]);
PD_p    = cumsum(p,3);
PD_p(PD_p==0) = NaN;